OpenStreetMap 結合 Python，如何找到你想要的資料¶
常用的 Python 套件與方法¶
slide: https://github.com/Supaplextw/COSCUP-2020/blob/master/Python-OSM-COSCUP-2020.ipynb
OpenStreetMap台灣 - 陳瑞霖
OpenStreetMap¶
- 地圖版的維基百科
- 在臺灣社群發展十年
- 自己的地圖自己畫
- tag 系統有彈性，能照顧少數需求
- 例：第三世界防災，台語地圖
運用套件和資料¶
- 直接用常用的網路套件：requests、json
- Overpass API
- 專用基礎套件：OSMPythonTools
- 進階套件：OSMnx
- 地圖呈現套件：folium
requests、json¶
- 簡單，只要有架設網路相關套件就能取用 OpenStreetMap
- 難度，需要有基本的網頁知識與 API 概念
OpenStreetMap 的 API¶
- OpenStreetMap API
- Overpass API
- 其他第三方 API，例如 Mapbox、Mapillary
In [1]:
import requests
import json
In [2]:
overpass_url = "http://overpass-api.de/api/interpreter"
overpass_query = """
[out:json];
area["nat_ref"="63"][admin_level=4];
(node["amenity"="restaurant"](area);
way["amenity"="restaurant"](area);
rel["amenity"="restaurant"](area);
);
out center;
"""
In [3]:
response = requests.get(overpass_url,
params={'data': overpass_query})
data = response.json()
In [4]:
import numpy as np
import matplotlib.pyplot as plt # 輸出圖表用
# Collect coords into list
coords = []
for element in data['elements']:
if element['type'] == 'node':
lon = element['lon']
lat = element['lat']
coords.append((lon, lat))
elif 'center' in element:
lon = element['center']['lon']
lat = element['center']['lat']
coords.append((lon, lat))
In [5]:
# Convert coordinates into numpy array
X = np.array(coords)
plt.plot(X[:, 0], X[:, 1], 'o')
plt.title('Restaurant in Taipei')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.axis('equal')
plt.show()
Overpass Turbo 範例：臺北市的餐廳 https://overpass-turbo.eu/s/WDS
匯出存檔為 geojson 檔案，直接取得 GeoPandas 能處理的格式
In [6]:
import geopandas as gpd
gdf = gpd.read_file("restaurant-tpe.geojson")
gdf.head()
Out[6]:
In [7]:
# Select some useful cols and print
cols = ['name', 'name:en', 'cuisine', 'addr:street', 'takeaway', 'delivery']
# Print only selected cols
gdf[cols].head(10)
Out[7]:
In [8]:
import osmnx as ox
import matplotlib.pyplot as plt
%matplotlib inline
ox.config(log_console=True, use_cache=True)
ox.__version__
Out[8]:
In [9]:
place = "Daan District, Taipei City, Taiwan"
G = ox.graph_from_place(place)
In [10]:
ox.stats.basic_stats(G, area=None, clean_intersects=False, tolerance=15, circuity_dist='gc')
# 傳回基本的統計資訊
Out[10]:
In [11]:
fig, ax = ox.plot_graph(G, fig_height=12, node_size=0, edge_linewidth=0.5)
# 傳回大安區範圍內的道路網
In [12]:
area = ox.gdf_from_place(place)
# 轉換 place 變成 GeoPandas 能處理的類型
In [13]:
type(area)
# 發現 area 是屬於 GeoPandas 的 GeoDataFrame 類型物件
Out[13]:
In [14]:
area
# 傳回大安區的最東、西、南、北的邊界
Out[14]:
In [15]:
area.plot()
# 用 Matplotlib 畫出臺北市的範圍
Out[15]:
In [16]:
area['geometry']
# 傳回臺北市的 geometry 數值，以及前幾筆構成 polygram 的點經緯度
Out[16]:
In [17]:
buildings = ox.footprints_from_place(place)
In [18]:
len(buildings)
# 回傳建築物數量
Out[18]:
In [19]:
buildings.head(3)
# 傳回前三筆建築資料，先大略看資料的情形
Out[19]:
In [20]:
buildings.columns
# 觀察 OpenStreetMap 資料當中，建築物有那些欄位
Out[20]:
In [21]:
# 截取餐聽
restaurants = ox.pois_from_place(place, amenities=['restaurant'])
# 大安區有幾間餐廳?
len(restaurants)
Out[21]:
In [22]:
# 有那些欄位
restaurants.columns
Out[22]:
In [23]:
# 選擇你感興趣特定的欄位
cols = ['name', 'name:en', 'cuisine', 'addr:street', 'takeaway', 'delivery']
# 只顯示特定欄位
restaurants[cols].head(10)
Out[23]:
osmnx¶
- 強項為路網分析
- networkx 衍生的套件
In [24]:
# 截取道路網的線條和交叉點
nodes, edges = ox.graph_to_gdfs(G)
In [25]:
nodes.head()
# 路網當中的交叉點資訊
Out[25]:
In [26]:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(18,12))
# Plot the footprint
area.plot(ax=ax, facecolor='black')
# Plot street edges
edges.plot(ax=ax, linewidth=0.6, edgecolor='#BC8F8F')
# Plot buildings
buildings.plot(ax=ax, facecolor='khaki', alpha=0.9)
# Plot restaurants
restaurants.plot(ax=ax, color='green', alpha=0.9, markersize=4)
plt.tight_layout()
In [27]:
from pyproj import CRS
# Set projection
projection = CRS.from_epsg(3826) #TWD97 投影
# Re-project layers
area2 = area.to_crs(projection)
edges2 = edges.to_crs(projection)
buildings2 = buildings.to_crs(projection)
restaurants2 = restaurants.to_crs(projection)
# 轉換座標系統，變成 TWD97
In [28]:
fig, ax = plt.subplots(figsize=(18,12))
# Plot the footprint
area2.plot(ax=ax, facecolor='black')
# Plot buildings
buildings2.plot(ax=ax, facecolor='pink', alpha=0.7)
Out[28]:
In [29]:
fig, ax = plt.subplots(figsize=(18,12))
# Plot the footprint
area2.plot(ax=ax, facecolor='black')
# Plot street edges
edges2.plot(ax=ax, linewidth=0.6, edgecolor='dimgray')
# Plot buildings
buildings2.plot(ax=ax, facecolor='pink', alpha=0.7)
# Plot restaurants
restaurants2.plot(ax=ax, color='yellow', alpha=0.7, markersize=4)
plt.tight_layout()
Wikidata 上的 OpenStreetMap 資料¶
- 幫忙打 Wikidata 廣告
- 介紹 Wikidata Query Service
- 也可以用 Python 來截取資料
In [31]:
# pip install sparqlwrapper
# https://rdflib.github.io/sparqlwrapper/
import sys
from SPARQLWrapper import SPARQLWrapper, JSON
endpoint_url = "https://query.wikidata.org/sparql"
query = """SELECT DISTINCT ?item ?name ?Hu ?OpenStreetMap ?Coor WHERE {
?item wdt:P31 wd:Q705296;
wdt:P131 wd:Q1867.
SERVICE wikibase:label {
bd:serviceParam wikibase:language "[AUTO_LANGUAGE],zh".
?item rdfs:label ?name.
}
OPTIONAL { ?item wdt:P5020 ?Hu. }
OPTIONAL { ?item wdt:P402 ?OpenStreetMap. }
OPTIONAL { ?item wdt:P625 ?Coor. }
}
ORDER BY (?itemLabel)"""
def get_results(endpoint_url, query):
user_agent = "WDQS-example Python/%s.%s" % (sys.version_info[0], sys.version_info[1])
# TODO adjust user agent; see https://w.wiki/CX6
sparql = SPARQLWrapper(endpoint_url, agent=user_agent)
sparql.setQuery(query)
sparql.setReturnFormat(JSON)
return sparql.query().convert()
results = get_results(endpoint_url, query)
for result in results["results"]["bindings"]:
print(result)
Python 的地圖套件：
- IPython.display
- Folium
In [34]:
%matplotlib inline
import ogr
import shapely
from shapely.geometry import *
import geopandas as gpd
import matplotlib.pyplot as plt
import IPython.display as display
osm = """
<iframe width="500" height="400" frameborder="0" scrolling="no" marginheight="0" marginwidth="0" \
src="http://www.openstreetmap.org/export/embed.html?bbox=121.5391%2C25.0175%2C121.5489%2C25.0106&layer=mapnik" \
style="border: 1px solid black"></iframe><br/><small><a \
href="http://www.openstreetmap.org/#map=17/25.01392/121.54314">View Larger Map</a></small>
"""
display.HTML(osm)
Out[34]:
In [35]:
import folium
m = folium.Map(location=[25.01508, 121.54275], zoom_start = 17)
m
Out[35]:
In [36]:
folium.Marker(location=[25.01508, 121.54275], popup="COSCUP 2020",tooltip="Click").add_to(m)
m
Out[36]:
In [37]:
folium.Marker(location=[25.01508, 121.54275], popup="研揚大樓",tooltip="Click").add_to(m)
m
Out[37]:
In [46]:
m.save('index.html')
參考資料¶
Reveal.js
Overpass API
OSMnx:
- OSMnx Example：路網分析
- 赫爾辛基大學 Geo-Python-AutoGIS 課程-Retrieving OpenStreetMap data
- OSMnx 說明文件
- Day20 OSMnx應用
IPython:
- http://darribas.org/gds15/content/labs/lab_03.html
- https://automating-gis-processes.github.io/site/notebooks/L5/interactive-map-folium.html
Folium:
- https://www.youtube.com/watch?v=QpBmO35pmVE
- https://github.com/tsaisunghao/Geography/blob/master/EarthquakeMap/geomapping.ipynb
- https://github.com/victorgau/wzu20180122/blob/master/%E5%9C%B0%E5%9C%96%E7%AF%84%E4%BE%8B/%E7%94%A8%20folium%20%E7%95%AB%E5%9C%B0%E5%9C%96.ipynb
- https://github.com/gotoariel/folium-demo/blob/master/demo/demo.ipynb
In [ ]: